﻿//--------------------------------------------------------------------------------------------------------------------------------
// <copyright file="SlaeGenerator.cs" company="Exul">
//     Copyright (c) Exul. All rights reserved.
// </copyright>
// <author>Exul</author>
//--------------------------------------------------------------------------------------------------------------------------------

namespace ExulLibrary.FiniteElementMethod.Statements.AxiallySymmetricRectangular2DScalarMaxwell
{
    using System;

    using ExulLibrary.Mathematics.LinearAlgebra;

    /// <summary>Provides methods to generate portrait of system of linear algebraic equations for the task.</summary>
    public sealed class SlaeGenerator
    {
        /// <summary>Generates the portrait of system.</summary>
        /// <param name="horizontalPointsCount">The count of points by the first dimension.</param>
        /// <param name="verticalPointsCount">The count of points by the second dimension.</param>
        /// <returns>The matrix of system with initialized portrait.</returns>
        /// <exception cref="ArgumentOutOfRangeException">
        ///     <paramref name="horizontalPointsCount"/> or <paramref name="verticalPointsCount"/> is less than 2.
        /// </exception>
        public CrcsSparseSymmetricMatrix<double> GeneratePortrait(int horizontalPointsCount, int verticalPointsCount)
        {
            #if SAFE

            if (horizontalPointsCount < 2)
            {
                throw FiniteElementMethod.ThrowHelper.PointsCountIsLessThanRequired(
                    "horizontalPointsCount", horizontalPointsCount, 2);
            }

            if (verticalPointsCount < 2)
            {
                throw FiniteElementMethod.ThrowHelper.PointsCountIsLessThanRequired(
                    "verticalPointsCount", verticalPointsCount, 2);
            }

            #endif

            var pointsCount = horizontalPointsCount * verticalPointsCount;

            // For rectangular grid an internal node has 4 relationships:
            // Y4  X  0
            // Y1 Y2 Y3
            var elementsCount = 4 * pointsCount;

            // nodes at left border don't have Y4
            elementsCount -= verticalPointsCount;

            // nodes at bottom border don't have Y2
            elementsCount -= horizontalPointsCount;

            // nodes at left and bottom borders don't have Y1 (-1 because left bottom corner belongs to both borders)
            elementsCount -= horizontalPointsCount + verticalPointsCount - 1;

            // nodes at right and bottom borders don't have Y3 (-1 because right bottom corner belongs to both borders)
            elementsCount -= horizontalPointsCount + verticalPointsCount - 1;

            var begins = new int[pointsCount + 1];
            var indexes = new int[elementsCount];

            // start index
            begins[0] = 0;

            var nodeIndex = 0;

            // left bottom corner node doesn't have relationships
            begins[1] = 0;
            nodeIndex++;

            // nodes on bottom border
            for (var index = 1; index < horizontalPointsCount; index++)
            {
                // nodes on bottom border have only Y4 relationship
                begins[nodeIndex + 1] = begins[nodeIndex] + 1;

                // relationship with Y4
                indexes[begins[nodeIndex]] = nodeIndex - 1;

                nodeIndex++;
            }

            // other nodes
            for (var verticalIndex = 1; verticalIndex < verticalPointsCount; verticalIndex++)
            {
                // node on left border have only Y2 and Y3 relationships
                begins[nodeIndex + 1] = begins[nodeIndex] + 2;

                // relationship with Y2
                indexes[begins[nodeIndex]] = nodeIndex - horizontalPointsCount;

                // relationship with Y3
                indexes[begins[nodeIndex] + 1] = nodeIndex - horizontalPointsCount + 1;

                nodeIndex++;

                // internal nodes
                for (var horizontalIndex = 1; horizontalIndex < horizontalPointsCount - 1; horizontalIndex++)
                {
                    // internal nodes have all four relationships
                    begins[nodeIndex + 1] = begins[nodeIndex] + 4;

                    // relationship with Y1
                    indexes[begins[nodeIndex]] = nodeIndex - horizontalPointsCount - 1;

                    // relationship with Y2
                    indexes[begins[nodeIndex] + 1] = nodeIndex - horizontalPointsCount;

                    // relationship with Y3
                    indexes[begins[nodeIndex] + 2] = nodeIndex - horizontalPointsCount + 1;

                    // relationship with Y4
                    indexes[begins[nodeIndex] + 3] = nodeIndex - 1;

                    nodeIndex++;
                }

                // node on right border have only Y1, Y2 and Y4 relationships
                begins[nodeIndex + 1] = begins[nodeIndex] + 3;

                // relationship with Y1
                indexes[begins[nodeIndex]] = nodeIndex - horizontalPointsCount - 1;

                // relationship with Y2
                indexes[begins[nodeIndex] + 1] = nodeIndex - horizontalPointsCount;

                // relationship with Y4
                indexes[begins[nodeIndex] + 2] = nodeIndex - 1;

                nodeIndex++;
            }

            return new CrcsSparseSymmetricMatrix<double>(begins, indexes);
        }
    }
}